# This is an detailed implementation of phylogenetic mark correlation functions.
# 
# data.ppp: a marked point pattern
# marks: list("mark"=?,"biva"=list("mark"=?,"from"=?,"to"=?))
# phyd: a symatric matrix of phylogenetic similarity distance, with
#       species rowname.
# exclude_sametype: wether exclude the samle type of points (e.g. the
#       conspecific individuals) in calculation.
# to:   maximum spatial distance to be calculated
# step: width of the cout interval
# scale: whether normalize the mark correlation funciton
#
# Author: guochun
###############################################################################



phylo_corr=function(data.ppp,marks,phyd,exclude_sametype=TRUE,
		to,step,neach,scale,report=TRUE){
	#remove NA
	sp=get_marks(data.ppp,marks=marks$mark)
	phyrow=rownames(phyd)
	spmatch=match(sp,phyrow)
	del=which(is.na(spmatch))
	if(length(del)!=0){
		data.ppp=data.ppp[-del]
		sp=sp[-del]
		spmatch=spmatch[-del]
		phyd=phyd[-del,]
	}
	
	#tranfer types of sp
	S=length(phyrow)
	sp=spmatch-1
	if(is.null(marks$biva)){
		to_points=data.ppp
		from_points=data.ppp
		fsp=sp
	}else{
		biva_mark=get_marks(data.ppp,marks$biva$mark)
		
		from_points=data.ppp[biva_mark==marks$biva$from]
		if(marks$biva$to!="all"){
			to_points=data.ppp[biva_mark==marks$biva$to]
		}else{
			to_points=data.ppp
		}
		
		fsp=sp[biva_mark==marks$biva$from]
	}
	
	#estimate normalized constant c
	if(scale)
		c=estimateC(sp,fsp,phyd,exclude_sametype)
	
	phyd=as.numeric(phyd);
	nloop=ceiling(from_points$n/neach)
	if(report)
		print(nloop)
	
	re=auto_lapply(X=1:nloop,FUN=mkc.core,fx0=from_points$x,fy0=from_points$y,
			fsp0=fsp,x0=to_points$x,y0=to_points$y,sp0=sp,phyd0=phyd,
			rmax0=to,step0=step,neach=neach,n=from_points$n,S=S,report=report)
	
	for (i in 1:length(re)){
		temp=re[[i]]
		if(i==1){
			totalbincout=temp[[1]]
			totalbinsum=temp[[2]]
		}else{
			totalbincout=totalbincout+temp[[1]]
			totalbinsum=totalbinsum+temp[[2]]
		}
	}
	
	result=totalbinsum/totalbincout
	
	if(scale){
		result=result/c
	}
	result=data.frame(r=seq(0,to,length.out=to/step),pk=result)
	attr(result,"constant")=c
	return(result)
}


mkc.core=function(i,fx0,fy0,fsp0,x0,y0,sp0,phyd0,rmax0,step0,neach,n,S,report=TRUE){
	if(report)
		print(i)
	if((neach*i)<=n)
		select=((i-1)*neach+1):(neach*i)
	else{
		select=((i-1)*neach+1):n
		neach=length(select)
	}
	lfx0=fx0[select]
	lfy0=fy0[select]
	lfsp0=sp0[select]
	
	nfocal=length(lfx0)
	ntotal=length(x0)
	bincout=rep(0,times=rmax0/step0)
	binsum=rep(0,times=rmax0/step0)
	
	temp=markcorr_innerC(lfx0,lfy0,lfsp0,x0,y0,sp0,rmax0,phyd0,nfocal,ntotal,step0,S,
			bincout,binsum)
	
	return(list(temp$bincout,temp$binsum))
}

markcorr_innerC <- cfunction(signature(fx="numeric",fy="numeric",fsp="numeric",
				x="numeric",y="numeric",sp="numeric",rmax="double",phyd="numeric",
				nfocal="integer",ntotal="integer",step="double",nsp="integer",
				bincout="numeric",binsum="numeric"),
		body=paste(readLines("./mark correlation function/markcorr_innerC.c"),collapse = "\n"),convention=".C",
		includes="#include <math.h>")


estimateC=function(sp,fsp,phyd,exclude_sametype){
	fab=table(fsp)
	tab=table(sp)
	if(exclude_sametype){
		diag(phyd)=NA
	}
	phyd=phyd[as.numeric(names(fab))+1,]
	phyd=phyd[,as.numeric(names(tab))+1]
	w=as.numeric(fab)*rep(as.numeric(tab),each=length(fab))
	phyd=phyd*w
	dim(w)=dim(phyd)
	diag(w)=0
	phyd=as.numeric(phyd)
	c=sum(phyd,na.rm=T)/sum(w)
	return(c)
}
